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Abstract 

We consider simple examples illustrating some new features of the linear response 
theory developed by Ruelle for dissipative and chaotic systems [J. of Stat. Phys. 
95 (1999) 393]. In this theory the concepts of linear response, susceptibility and 
resonance, which are familiar to physicists, have been revisited due to the dynamical 
contraction of the whole phase space onto attractors. In particular the standard 
framework of the "fluctuation-dissipation" theorem breaks down and new resonances 
can show up oustside the powerspectrum. In previous papers we proposed and 
used new numerical methods to demonstrate the presence of the new resonances 
predicted by Ruelle in a model of chaotic neural network. In this article we deal 
with simpler models which can be worked out analytically in order to gain more 
insights into the genesis of the "stable" resonances and their consequences on the 
linear response of the system. We consider a class of 2-dimensional time-discrete 
maps describing simple rotator models with a contracting radial dynamics onto the 
unit circle and a chaotic angular dynamics Ot+i = 20f(mod 27r). A generalisation 
of this system to a network of interconnected rotators is also analysed and related 
with our previous studies [8,11]. These models permit us to classify the different 
types of resonances in the susceptibility and to discuss in particular the relation 
between the relaxation time of the system to equilibrium with the mixing time 
given by the decay of the correlation functions. Also it enables one to propose some 
general mechanisms responsible for the creation of stable resonances with arbitrary 
frequencies, widths, and dependency on the pair of perturbed /observed variables. 
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1 Introduction 



In statistical physics, a standard idea is that the hnear response of a system can be com- 
puted in terms of correlation functions of this system at equilibrium. This is the cornerstone 
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of important results like the fluctuation-dissipation theorem [1,2]. More recently a general 
theory for the linear response for hyperbolic systems has been developed by Ruelle [3,4]. 
This theory shows that the picture of correlation functions is incomplete in this case. In this 
situation a linear response function has in general two contributions, corresponding respec- 
tively to the expanding (unstable) and the contracting (stable) directions of the hyperbolic 
dynamics. Although the first contribution, called in the sequel the unstable response func- 
tion, can still be associated with some correlation function, this is not true for the second 
contribution, named the stable response functiocQ- This property transfers to the complex 
susceptibility, i.e. the frequency-dependent responses of the system, for which we can also 
distinguish between two types of resonances, namely the stable and the unstable resonances. 

Computations of the linear response in specific chaotic systems have been performed by 
several authors [5]- [7]. For example in [5], Falcioni et al. discuss and numerically apply a 
simple method to calculate the (impulse) linear response function. The procedure is an 
averaging on perturbations of the system by means of Dirac kicks. However, as estimated 
by the authors, the convergence of this method is slow in the time-domain. Another study 
has been performed by Reick in [7], who uses a different technique, averaging the effects of 
time-periodic perturbations. With this method he numerically demonstrates the existence 
of linear response in the Lorentz system, although the latter is known to be nonhyperbolic. 
Neither of these previous works, however, obtain evidence of the stable resonances predicted 
by Ruelle. Let us remark that contrary to the present study, these papers do not consider 
chaotic systems with state-dependent perturbations. 

In a recent work we have illustrated for the first time the prediction of Ruelle on a simple 
model of networks of interconnected units exhibiting chaotic dynamics [8]. Using indepen- 
dently the same averaging procedure as [7], we managed to compute the susceptibility Xijiyj) 
and its complex poles. (Here Xiji}^) denotes the complex amplitude of the average response 
of unit i to a harmonic excitation of unit j at frequency uj). Our numerical method could 
not provide separately the unstable and the stable contributions in the susceptibility, but 
our results demonstrated the existence of stable poles in Xijif^) by setting in evidence poles 
which were not in the set of poles of the correlation functions (the Ruelle-Pollicott reso- 
nances [9,10]). 

On the other hand our numerical results pointed out towards interesting features whose 
consequences are worth to explore further. First, stable resonances in Xij can be specific to 
the pair ij. This property, which is not true in general for the Ruelle-Pollicott resonances, 
is quite useful for networks as it can be used to transmit an amplitude modulated signal 
between two specific units in the network [11]. Another unexpected feature is that the average 
relaxation time of unit i to an impulse perturbation of unit j (which can be predicted by 



Though the terminology "stable", "unstable" is a little bit confusing when dealing with hnear 
response and resonances, we use it because it refers explicitely to unstable and stable spaces in the 
tangent space of the attractor. The mechanism leading to the existence of a stable and unstable 
contribution of the linear response are of completely different nature: this is exponential contraction 
for the stable part and this is exponential mixing for the unstable part. 
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Fourier transform of Xij{uj)) can extend on a much longer time lapse than the mixing time, 
i.e. the decay time associated with the correlation functions. 

These properties cannot be understood in the framework of the classical fluctuation-response 
theorem. Indeed they reflect the presence of the stable contribution of the response functions 
which arc not correlation functions. However the stable response or the stable susceptibility 
are difficult to extract numerically. Therefore the goal of the present paper is to study a class 
of toy models for which the stable and unstable decomposition of response functions can be 
worked out explicitly. In particular these examples clearly show the relative contributions of 
correlation and non-correlation terms in the response functions, and their treatment helps 
to better understand the above properties revealed in networks. These simple examples are 
also treated numerically. This allows one to validate the numerical procedure used in [8,11]. 

In the next Section we recall some elements of Ruelle's theory which will be extensively used 
later on. We formally derive his general response formula by using the method of the impulse 

response and then introduce the definition of susceptibility as being its Fourier transform. 
In Section 3 we treat the basic example 6t^i = 26t (mod 27r) in this framework and discuss 
briefiy the issue of the non-invertibility of this dynamics. This elementary example serves also 
to describe and to test numerical methods. Section 4 concerns class of uniformly hyperbolic 
systems which modelise chaotic rotators in the plane. Various extensions of the simplest 
model are considered in the following subsections so as to discuss properties of the stable 
resonances. They are summarised in the Conclusions. An appendix of the paper presents a 
natural extension of our chaotic rotator model in higher dimensions for which some properties 
generalise directly. 



2 A general linecir response formula of Ruelle 

We start by reviewing some results of the linear response theory obtained by Ruelle [3,4]. We 
concentrate on the particular case of autonomous dynamical systems described by iterations 
of a map. Consider the following systems whose dynamics is governed by the recurrence 
equation: 



where t G Z, x G M, and M is the phase space, e.g. a compact manifold in M . The map 
is defined by a smooth function F(x), not necessarily invertible, and Xj denotes the state of 
system at time t. We also use the notation F(xt) = F*(xo). The dynamics of (1) is assumed 
to be chaotic, mixing and associated with an ergodic measure pp of Sinai-Ruelle-Bowen type 
(SRB) [4,12]. This implies that any measurable observable ^(x) has a time average 



Xt+l 





_ 1 ^ 

lim -y A(F\x)) 



(2) 
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which is equal to its ensemble average: 

<A>= J pj.(dx)yl(x) = ^lim J dx A(F*(x)) (3) 

for almost every initial conditions x selected with respect to the Lebesgue measure (ix. Note 
that in eq. (3), pp equivalently denotes the image of the Lebesgue measure under F*, i.e. 
Pf = limj^oo F*^ diL where the system is assumed to be mixing. 

Example: 

We consider a first model which will be dealt with in next Sections. This simplest example 
for system (1) is given by the mapping: 

^t+i = f{Ot) (4) 

with f : 6 ^ 26 (mod 27r) defined on the unit circle. This system is chaotic with Lyapunov 
exponent log 2. This means that a small perturbation 69 of any initial condition 6 is locally 
amplified in time with speed 2*. On the other hand, on the unit circle, the difference between 
the initial trajectory and the perturbed one is blurred in time. In other words, despite the 
sensitivity to initial condition, the mean effect of the perturbation 59 applied at time t — 
only on observables of this system vanishes in time. This is readily seen as a consequence of 
the ergodic property A —< A >, namely: 



for almost every initial conditions. Let us recall that in this example the SRB measure is 
d9/2n and so < A >^ ^ J^^'dOAie). 

Therefore in the class of chaotic systems considered here, a single perturbation of initial 
condition has no effect on the mean value of any observable, although the short time effect 
of this perturbation can be important. Now let us suppose that one applies on system (1) 
a permanent perturbation depending on time t and possibly on the state of this system. 
Then the goal of the response theory is to study how the mean values of observables are 
affected by this change, compared with the unperturbed system. This problem is not easy 
in general because in this case the SRB state is asymptotically different from the one of the 
unperturbed system. 

Thus now let us consider a perturbed version of system (1) governed by the following equa- 
tion: 

i,+i=F,(i,)4^fF(i,)+^,+,(x,) (5) 
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The term (x) defines the perturbatioiH]- The 
at time s will be denoted Xf = F^^'^^Stg)- (In the 
on X when no confusion is possible). Then one 
with respect to the perturbed system as [3]: 

< A>t= ^hrn^J rfxA(F(*'')(x)) 



iterate at time t of initial condition x^ taken 
following we will abandon the tilde notation 
can define the mean value of observable A 

(6) 



2. 1 The impulse linear response 



The goal of the linear response theory aims to compute < A >t — < A > to first order in ^. 
This will be denoted by < 5 A >t. The relevance of this concept in chaotic dynamical systems 
has been analysed by Ruelle. The latter introduced the concept of differentiation of a SRB 
state in [13] and provided general formula for its computation in [3,4] under the asumption 
of uniform hyperbolicity. 



We (formally) derive the formula of Ruelle in a different way by using the method of the 
impulse perturbation. Consider a perturbation of the form ^^(x) = ^(x) Str- In this case, and 
if t > r > s one can write: 



= o . . . o . . . o ^ pt-r o (F + ^) O F- 



r— s — 1 



since Fj = F for all t except for Ft—i = F + ^. Therefore eq. (6) becomes: 



<A>t=\im I rfx \a o F*-^ o (F + ^) o F 

s^-co J L 



T — S — 1 



AoF 



t-T 



[F(x)+^(x)) 



where we have used eq. (3) to take the limit s — oo. On the other hand, by using the 
invariance of pp under F, one can write < A >= J pp^d^s.) [A o F*^^] (F(x)) . Consequently, 
the impulse response at time t > t can be set under the following form: 



< A>t- < A>= / p^(rfx) {A o F*-")(F(x) + ^(x)) - (A o F*-0(F(x)) 



(7) 



Although the expression (7) is valid beyond the linear regime, this form is useful for com- 
puting the linear response < 6 A >t. The latter will be denoted XAi{t — t) and can be readily 

^ The subscript t + 1 in eq. (5) is neither a missprint neither a violation of causality. The motivation 
to denote the perturbation acting on xt by instead of is that it allows one to write the 
subsequent convolution formula (9) and the Fourier transform (10) into the standard forms they 
have for continuous time systems 
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deduced as the first order Taylor expansion of (7): 
XA^t - ^) = / Pf(c^x)V(^ o F*--)(F(x)) • ^(x) 



(8) 



Now a general time-dependent perturbation £j(x) can be written as the sum of impulses 
^t(x) = Z]r=-oo ^t(^) <^tr- Thus, the linear response is given by a superposition of terms 
like (8), namely: 

<5A>,^ ^ XAi{t-r)-U^) (9) 

T= — CX> 

This general convolution formula is the linear response derived by Ruelle in [3,4]. Note 
that it is not evident that the series converge. This point was initially raised in the famous 
objection of van Kampen [14]; Due to the existence of unstable directions and positive 
Lyapunov exponents, one would actually expect that the linear response diverges. However, 
this issue was highlighted by the statistical physics arguments of Kubo [15] and by several 
subsequent studies, as the ones of Falcioni et al. [5] and the recent contribution of Bofetta et 
al. [6]. In the present context where one assumes the uniform hyperbolicity, one can follow 
the argument of Ruelle who has shown, by using projections on the local stable and unstable 
subspace, that the unstable part is a correlation function, thus the series converges due to 
(exponential) mixing occuring in uniformly hyperbolic systems; on the other hand the stable 
part converges due to exponential contraction. 



2.2 The susceptibility 

From eqs. (8)-(9), considering perturbations of the form ^^(x) = ^(x) e~*'^*, it is easy to show 
that the linear response of (the mean value of) observable A is given by 

<5A XAd^)e-'^' 
where the complex amplitude 

oo 

XmH = Y: XMit) e'^' (10) 

t=—oo 

is called the susceptibility. The latter is nothing but the Fourier transform of the impulse 
response XA£,{t) of observable A with respect to perturbation ^ (with XA£,{t) = for t < ). 
Conversely, the impulse response can be computed from the susceptibilty by inverse Fourier 
transform XAi{t) = ^ Jq'" XAi{u^) e'^^du. 
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Our derivation of eq. (9) gives a straightforward but non rigorous way to obtain the hnear 
response, and in particular the susceptibihty, for dynamical systems of the form (1). The 
rigorous analysis has been done by Ruelle. In fact, in refs. [3,4] using the hypothesis of 
uniform hyperbolicity of the systems. Ruelle goes much beyond the eqs. (8)-(10) by proving 
a general decomposition of the response functions into two terms -stable and unstable- 
which stems from the natural foliation of the phase space into stable and unstable manifolds 
and its implications for the SRB measure. In the following Sections, instead of rephrasing 
Ruelle's theory about this decomposition proved for a general system, the goal is to apply this 
concept on simple systems in order to figure out some of its consequences already observed 
in our previous works [8,11]. Prior to analyse these effects, we review in next paragraph the 
simplified case where the system has only expanding directions and so no contractions in 
phase space. 



3 Purely expanding dynamics. 

3.1 A basic example. 

We start by considering the model of example 1 introduced above {f : 6 ^-^ 2^ mod 2tt), 
perturbed by a small periodic forcing e^t = e^e~*'^*, namely: 

9t+i = fiet) + ea0t)e-'^^'^'^ (11) 

Let us consider the observable A{9) = "the value of 9 in [0, 27r)"[l]. When e = 0, we have 
< A{9) >= TT . Then, following the discussion of the previous Section, the mean value of 
A{9) in the perturbed system is given by < A{9) >t= vr + ex{uj)e~'^'^^. In this expression the 
susceptibility x{uj) is the Fourier transform of the impulse response xit) ^-nd the latter can 
be computed according to eq. (8) as: 

x{t) = -jd9de{Aof%f\e))m (12) 



for t > 0. The notation d0{A o P){f{9)) means the derivative w.r.t. 9 of A{f^{9)) evaluated 
at f{9). To treat eq. (12) one wishes to use a change of variable = f{9) followed by an 
integration by part. This would lead us to the fluctuation response theorem. Since / is not 
assumed invertible, the proposed change of variable is not straightforward. We shall come 

^ Let us remark that whereas the observable A{9) is a natural choice to measure the angle, it is 
discontinuous for each = nlir (n integer), with a discontinuity equal to — 27r . Thus the derivative 
of A standing in eq. (8) should be treated as deA = 1 — 27r ^„ 6{6 — nln). In fact, as it is shown 
below, the method of integration by part avoid the explicit use of this distribution. 
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back to the general case later on. For simplicity we first consider a particular situation, where 
the perturbation function ^(6') can be written as {X o f)[6). In this case, a change of variable 
= f{0) can be performed in (12). This change of variable can be followed by an integration 
by part leading to the following expression (we return to the variable 6): 

Xit) = X(27r)5,,o -^JdO Aifie)) deXiO) (13) 





The first term in the right-hand side of this equation results from the integrated term 
^[A{f{9))X{e)]l'^. This is consistent with the trivial case where X is a constant. Then 
the impulse response is zero except at t = where it is equal to this constant. Conversely 
the integrated term disappears for an observable A and a perturbation X continuous in 6. In 
this case one obtains a clearcut application of the fluctuation-response theorem: the linear 
response is indeed a correlation functior0 (truncated to for t < 0). This can be expressed 
as: 

x{t) = -<9t;dgX{9o)> (15) 



for t > 0. As a concrete example , let us consider the function 

x{e) = j£y^e{2n-e) (i6) 



(The number (6/27r)^ is just a normalisation factor such that < X >= 1.) Then a simple 
calculation yields the following function: 

x{t) = 12<et;eo>=2-' (17) 



for t > 0. Indeed for model 1, all the correlation functions can be computed analytically [9]. 
In particular <9t;9o >= 2^712. 

Now we can also compute the susceptibility given by the Fourier transform of (17). This 
gives: 

, 1 2(2 — costu + i sino;) , , 

= i_ei(c.+iiog2) = 5_4cosa^ ^ ^ 



^ If j4(x) and -B(x) are two observables, the time correlation function of these observables is defined 
by: 

< A(xi); S(xo)) >= / pF(dx)A(F*(x))S(x)- <A><B> (14) 
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Let us remark that xi^) has a complex pole in cu — —i log 2 on the imaginary axis. 

3.2 Numerical methods 

The two elementary analytical results (17)-(18) can be used to test numerical methods to 
compute x{t) or x(c<j). In particular it can serve to validate the procedure we have used in 
analysing the linear response of a chaotic network of interconnected units [8]. 

Let us first consider a direct method to compute the impulse response introduced 
by Falcioni et al in [5]. This method was proposed for computing the impulse response to 
perturbations independent of the state, but it is easily adapted here where the perturbation 
depends on the state. In what follows we discuss the method for the system 6*4+1 = fiOt) 
introduced above, but of course the method is general. As before the goal is to compute 
< 5 A >t to first order in e, when the perturbed system follows the same dynamics as the 
unperturbed one, except that at time t — Q the state is displaced by ^0 = ^0 + e-'^(^o)- 

The principle of the method studied in [5] can be described as a finite-sampling estimation of 
the integral defined by eq. (7). More precisely, for i > fixed, is estimated by averaging 
many reahsations of A(f*{9k)) {k — 1, - ■ • ,N), starting from initial data for 9^ chosen with 
the SRB measure on the (unperturbed) attractor. Adapted to the present notations, the 
authors of [5] propose to evaluated x(i) as follows: 

1 ^ 

xit) «^ T7- E + ^^(^^)]) - A (fm)) (19) 

In practise, reckoning with the ergodicity of one trajectory of the unperturbed dynamics, 
the set {9i,92, ■ ■ ■ ,0n} can be prepared by considering only one (long) trajectory < t <T 
and dividing T in N samples, T — Nt, where r is the maximal time for which x(i) will be 
caculated. Then the initial data 9k can be simply chosen as f'^{9k-i). 

Now, considering the concrete example for X{9) worked out in last Section, figure 1 compares 
the response function calculated numerically with this method, with the analytical result of 
eq. (17). It shows a good agreement of the numerics for relatively short times, beyond which 
error fluctuations increase dramatically. In ref. [5] the authors discuss the growth of the error 
due to flnite size sampling and show that it grows faster than exp(At)/V^ where A is the 
largest Lyapunov exponent (here A = log 2). This difficulty is enhanced in computing the 
susceptibilty x(a;), as this Fourier transform needs a long time estimate of (19). Therefore 
in [8] we have developed a new method which takes the other way round. We compute 
numerically the susceptibility and, if needed, deduce from it the impulse response. 

The principle of the new procedure we have implemented is quite simple [16]: returning to 
the system (11), we know that the average < A >t differs from < A >= Trby < SA >t— 
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5 10 15 20 



t 

Fig. 1. The impulse response of system (11)-(16) is given by x{t) = 2~* (solid line). It is compared 
with numerical simulations using the estimate (19) (with e = 0.01, = 4552000, symbol "*" ). In 
this example, it is seen that large fluctuations grow from t = 8. The other points ( symbol "x" ) are 
obtained by using the indirect method which consists in calculating the (inverse) Fourier transform 
of estimate (21) (See also Fig. 2). 

ex(c^)e~*'^*. So x^uj) =< 6 A >t e*'^*/e is independent of t and we can write for arbitrary large 
T: 

xM = ;^E<M>,e-* (20) 
^ t=i 



Let us remark that in this equation < 5A >t=< A >t — < A > represents a complex 
number because the perturbed dynamics given by (11) is written by means of the complex 
perturbation e^(6')e~*'^*. In practise, however, i.e. for numerical purpose, < A >t can be 
decomposed in real and imaginary parts, < A >t=< Ai >t +i < A^ >t, and both of 
them can be computed separately by considering respectively perturbation of the forms 
eC,{0) cos{ujt) for < Ai >t and of the form —e^{9) smiut) for < A2 >t- Now the important 
point is that the above average < >t could be attained also by making use of an ergodic 
assumption. In other words, for a fixed uj the time-periodic SRB measure underlying the 
average < >t can be created by the long time evolution of the perturbed system. Therefore 
we can estimate (20) by writing, for 7^ 0, small e and sufficiently long T: 

xM^;^|:A(/*(e))e^-* (21) 
^ t=i 



where / is the perturbed dynamics given by (11), modulo its practical implementation in 
terms of real and imaginary parts. Let us notice that in the latter sum we can replace 5 A by 
A because the contribution of the non-perturbed dynamics in this average disappears (for 
u = 0). 
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Fig. 2. Modulus of the complex susceptibility as computed analytically [from eq. (18), solid line] 
and numerically [symbol "x"]. The latter is given by the estimate (21) averaged over 200 samples 
obtained with e = 0.01 and T = 2^^. 

Our numerical methods allows one, in principle, to compute the susceptibility x{uj) for an ar- 
bitrary sampling of u = k27T/M G [0, 27r]. The subsequent inverse Fourier transform can then 
provide the impulse response on the time interval [0, M] with large M. This procedure has 
been used with success in [8,11]. Here we illustrate it on the basic example considered above. 
Figure 2 shows comparison of the analytical result of eq. (18) with the numerically computed 
susceptibility, which is quite satisfactory. Note that both approximation relations (19)-(21) 
do not involve an explicit linearisation (or a differentiation) of the dynamics. Therefore 
these relations can be used to numerically tackle the nonlinear response of the system to 
perturbations. From the operational point of view, the linear regime is characterised by a 
non-dependence of the results in e. 

NB: To obtain these numerical results we had to iterate the dynamics 6t+i = 26'(mod 27r and 
it is well known that a naive simulation of this map eventually drives all trajectories to 0. 
We have circumvented the problem by encoding numbers with an (arbitrary long) sequence 
of bytes and used the left shift << available in C language. This shift set the last bit to 
zero. Then we replace it by a bit or 1 randomly chosen with a probability p or q. For 
p = q = ^ we obtain then trajectories typical for the SRB measure. Note that we can obtain 
trajectories typical for any Bernoulli measure by changing the probabilities p, q. We believe 
that this method is well known, though we haven't been able to find any reference dealing 
with it. 



3.3 Dealing with the non-invertibility of f 

Now we come back to the general case where the perturbation C,{0) cannot be written as 
{X o f){6). Then, the change of variable = f{6) used in the integral (12), can still be 
performed, provided that the integration domain [0, 2n] is decomposed into the subintervals 
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[0, tt] and [tt, 27r] where / is invertible. In this way the same formal expression (13) for 
is obtained as before, but where: 



x{e) 



(22) 



Thus in the chosen example, the effective perturbation X{9) appears as the perturbation ^ 
averaged over the pre-images of /. This concept generalises to general one-dimensional map 
/ defined on an interval I — [j^Ik assuming the restrictions of / to each invertible (say 
fk) by writing : 



where Xhi^) indicatrix of interval Ik] = C ° fe'' ^-^d Vk = V ° fe'^i assuming 

that the SRB measure is written Pf{d9) = r]{6)d9. Again a concrete example is easy to 
work out in the case of model 1. Let us choose ^{9) = 9{2'k — 9). Then eq. (22) becomes 
X{9) = 9{2t: -9)/2-l/A and it is found that x(^) exactly as eq. (17) but divided by 4. 



4 Chaotic rotator in the plane 



In the previous Section we have seen that in the case of systems with expanding dynam- 
ics only the fluctuation-response theorem applies: under suitable hypothesis, the response 
functions are matched with correlation functions. In view of the simple examples treated 
above, the technical point which makes the result work is integrating by part in eq. (8) 
(and in example eq. (12)). In the general case this integration by part is no longer possible 
(in the stable directions) for hyperbolic dynamics with contraction in phase space. This is 
because the SRB measure becomes usually singular in the stable directions so that it cannot 
be written with a density function, e.g. pp{d'x) ^ ri{x.)d:K for some regular function ?7(x). 
Nevertheless the response function given by the integral (8), and so the susceptibility, can 
be well defined, but in general cannot be identified anymore with a correlation function. 
To analyse this situation, Ruelle decomposes the perturbation giving rise to the response 
functions in two parts, by locally projecting this perturbation respectively onto the stable 
and the unstable directions of the tangent dynamics. This can be achieved when the system 
is uniformly hyperbolic, i.e. when the tangent space at any point x decomposes uniformly 
with respect to x into a direct sum of a stable and an unstable subspace (e.g. [9]). 

The goal of this Section is to illustrate this decomposition on simple models where the 
uniform hyperbolicity holds by construction and where explicit computations can be done. 
To this purpose we consider a class of two-dimensional hyperbolic systems with a chaotic 
attractor on the unit circle, as represented on figure 3. Here the phase space is represented by 
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the plane which is contracted by the dynamics onto the unit circle where there is a rotation 
with expanding dynamics. This might be one of the simplest example representing a chaotic 
(non-fractal) attractor. There are several ways to implement this qualitative dynamics, the 
simplest one being to decouple the rotation and the radial contraction. In all the forthcoming 
implementations of this class of model, one assumes that the SRB measure can be factorised 
as a product of the form: 

pp{dxdy) = pp{r-drd6) = S{r — l)r dr rj{6)d6 /2t{ 

y 

Jl cliaotic rotation 
raditd contracHon 



Fig. 3. Schematic representation of the phase portrait of the toy model. The phase plane is shrinked 
by the dynamics on the unit circle playing the role of the chaotic attractor. Models implementing 
this scheme are given by eqs. (28). 

Here 5{r—l) is the singular part of the measure, due to the contraction of the radial direction. 
On the other hand the measure on the unit circle (which represents the unstable manifold 
of any point on the attractor) is continuous and described by the density function ri{9). The 
latter is equal to 1 in the particular case of the dynamics 9 ^ 29 considered in the previous 
Section. Before specifying some implementations of the dynamics of Fig. 3, we summarize 
the objective of the linear response theory and introduce the notations as follows. Consider 
the 2-dimensional time-discrete systems: 



xt+i = F^{xt, Vt) + (-X ix{t + 1, Vt) 

yt+i = Fy(xt, yt) + ey ^y(t + 1, xt, yt) (23) 

where the first term F = {F^, Fy) is the (autonomous) dynamics which is perturbed by a 
time-dependent vector field of small amplitude (e = |(ea;,ey)| <^ 1). Here the coordinates 
[x, y) do not necessarily mean the cartesian coordinates, but general coordinates on the 
plane. In the following, wc will consider only two cases, however, namely the cartesian or the 
polar coordinates. Furthermore, the observables will be the coordinates themselves, except 
for the "angle" variable of the polar coordinates which will be chosen like previously as A{9), 
i.e. the projection of 9 in [0,27r). 
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Now suppose that in absence of perturbation the mean values of x and y are zero. Then, 
when the perturbation is turned on, we are interested to the mean values of {xt,yt) to first 
order in e. When ^(t,x,y) = ^(x, |/)e~*'^*, the solution of this problem has been discussed 
above and it can be written in terms of the susceptibility matrix, namely: 





fex] 




1 













(24) 



The matrix elements Xiji'^) {hj = x or y) are the Fourier transform of the impulse response 
matrix elements given by eq. (8) : 



Xijit) = J pFidxdy)\/F^iFix.)) ■ (e,(x)e,) 

PF{dxdy)VFl{^) ■ (X,(x)ej) (25) 



In the last equation the function Xj coincides with S,j o if the dynamics is invertible. If F 
is not invertible then, as discussed in the previous Section, Xj can be interpreted as a type 
of averagin^^] of over the pre-images of F. 

Finally, assuming that the system is uniformly hyperbolic, i.e. when the tangent space at 
X decomposes uniformly into a direct sum of a stable and an unstable subspace, the vector 
Xj(x)ej appearing in (25) will be decomposed as 

X,(x)e,=X^^ + XJ (26) 



Then, by substituting this sum in eq. (25) we will be mostly interested in the ensuing 
decomposition of Xijif) iiito the stable and unstable susceptibilities as follows: 



x.,w=xl;Ht)+xlrw (27) 

for t > (and Xij = otherwise). 

In what follows we start by considering the cases where by construction XJ = or X^ = 
in eq. (26). This enables one to concentrate on the stable or on the unstable susceptibilities 
and their resonances, the former being the novel feature of out of equilibrium dissipative 
systems. 



^ In the non-invertible case, a particular case where the computation of this averaging over the 
pre-images of F gives a well controlable result is when the perturbation function is chosen under 
the form = Xj oY. Then the result is simply Xj. 
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4-1 Stable resonances 



As a first realisation of tlie phase portrait depicted on Fig. 3, we use polar coordinates (i.e. 
X = 6 and y = r in the eqs. (23)) and we consider tlie following system: 



n+i = Rirt) + erUt + hOt,rt) (28) 

with g(9) = 29 + b{r - 1) (6 constant) and R{r) = 1 + e"^(r - 1), with /i > 0. 

We first rule out the case 6 = 0. Then the angular dynamics is the same as before, g{9) = 29, 
and the system becomes extremely simple because the radial and the phase dynamics are 
decoupled. The radial dynamics is trivially contracting on the unit circle, with = -R*(ro) = 
1 + e"^*(ro — 1). Also, the impulse response matrix is diagonal, because Xr6»(^^) = Xeri^^) = 
which is readily deduced from eq. (25). Using the same equation, the diagonal elements are 
computed as follows: 



Xee{t) = J PF{rdrd9) def{9) Xe{9, r) = -<9t; dgXg{9o, 1) > (29) 
Xrr{t) = J PF{rdrd9) drR\r) Xr{9, r) = e"^* < X,(l, 9) > (30) 

for t > 0. The first line is identical to the response function (15) obtained for the basic model 
studied in the preceding Section. This unstable response function, corresponding to perturba- 
tion parallel to the attractor, is thus a correlation function. As an illustrative example we can 
calculate these integrals in the case of perturbation function ^r(t, 9, r) = C,g(t, 9, r) = X{9)6to, 
where X{9) is chosen as before (cf. eq. (16)). Then the susceptibility xgg{uj) is given explicitly 
by eq. (18). One notices that it possesses a pol^ inu = —i log 2. On the other hand, eq. (30) 
is an example of stable response, i.e. the response to a perturbation along the contracting 
direction of the attractor. Note that the result is calculated without integration by part. The 
corresponding susceptibility is then computed as: 

C 

*-(^) = 1 _ (31) 



where C is a constant. Now this function has a pole located on the imaginary axis in u; = — i/x. 
This pole is independent of the ones of the correlation function. In particular it can be 
situated arbitrarily close to the real axis. This simple example clearly shows the effect of two 
types of perturbation on the system, either transverse or parallel to the attractor. 



" More generally, the Ruelle-Pollicott resonances of this map are given by — iA;log2 where A; is a 
positive integer [9]. 
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Let us deal with the case b ^ 0, which introduces a skew dynamics between the stable and the 
unstable variables. Then, if one starts with an initial condition (^O)^o) the angular dynamics 
is easily worked out, giving the following expression : 



0t = 9\0, r) = + b{r - 1)^3^) (mod 27r) (32) 

Let us notice that the dynamics on the attractor (r = 1) is the same as before. This implies 
that the diagonal terms of the susceptibility remains the same as for 6 = 0, and the correlation 
functions do not change. However, the dependence of on r creates now a non-diagonal term 
which can be written as follows: 

2* - p-'** '^r 2* 27r 

Xer{t) = b- —B{t) with B{t)= X{e)de-Y.X{k—)2-' (33) 



Let us remark that the factor B{t) can be interpreted as an "error" function measuring the 
difference between the integral of X{6) and its approximation by a Riemann sum with inter- 
vals 2~*. (This factor B{t) stems from the distribution derivative of the "angle" observable 
eq. (3) ). It is clear that if X[9) is Riemann integrable, then the "error B{t) has to decrease 
to zero when t becomes large. Actually, in view of eq. (33) it should decrease faster than 2~* 
to get Xer{t) going to with t large. For example in the case of X{9) given by eq. (16), we 
obtain B{t) — 2~^* and thus 

n-t _ p-(21og2+//)t 

Xerit) = b ^—^^ (34) 



This geometrical interpretation of function B{t) shows directly that if X{9) varies much 
in [0, 27r] then it will take a long time t to reach an actual exponential decrease of B(t). 
By the meantime the response function can transiently grow a lot due to the exponential 
factor 2*. Therefore the degree of variation of X{6) determines the transient growth of the 
stable response function. This property will be used later in Section 4.4 by considering a 
perturbation function of the form X{9) — cos(2^^). 

Now, the susceptibility xori^), which can be readily deduced from (34), possesses a pole 
situated in a; = —i (2 log 2 -|- n) on the imaginary axis. The /j, contribution is clearly due to 
the contracting dynamics, whereas the first term 2 log 2 is in fact a Ruelle-Policott resonance 
(in the present case they are multiples of log 2). Thus this example shows a new feature 
which can appears when the variables are not trivially decoupled into stable and unstable 
variables: some stable resonances become a (linear) combination of poles due to the con- 
tracting direction and of resonances belonging to the correlation functions. An attempt of 
generalisation of this property is given in the Appendix. 
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4-2 Stable resonances with arbitrary frequencies 



The elementary model described by eq. (28) has only a one-dimensional radial dynamics. 
Consequently any stable pole of the corresponding susceptibility will be located along the 
imaginary axis. On the other hand it is possible to generalise this model in order to show 
that stable poles can have arbitrary frequencies and relaxation times. It suffices to increase 
to number of degrees of freedom of the system. To illustrate this point we present in this 
section two examples for which the stable resonances can have arbitrary non zero frequencies. 

A simple way to add more degrees of freedom to the radial dynamics is to transform its 
dynamics into a higher order reccurence. For instance let us consider a new type of rotator 
whose equations are a variation of system (28) taking the form: 



Ot+i=g{Ot,rt) 

rt+i = R{rt,rt-i) (35) 

with as before g{e) = 26 + b{r - 1) but now R(r, g) = 1 + 2e"''cos(a;o)(r - 1) - e~'^'^(Q - 1) 
(;U > 0). The choice of this dynamics is motivated by the fact that it creates in the stable 
susceptibility a pole placed simply at ujq — ifi. Thus we can easily control its frequency by 
varying ujq. This example is easy to deal with because the radial dynamics can be solved 
exactly, and it can be checked that given tq arbitrary and r_i = 0, one gets: 

rt = 1 + (ro - l)e"''*sina;o(t + l)/sina;o (36) 



This case will be considered below in cartesian coordinates. 

Another way to increase the number of degrees of freedom is to consider a network of inter- 
connected rotators. This example is interesting because it has links with the chaotic network 
model that we have studied previously [8,11] and which is mentionned in the introduction 
of this paper. So we consider now a collection of N interconnected rotators of the same type 
as before (cf. fig. 3), assuming for simplicity that the coupling between the units occurs only 
by means of their radial variables. The equations of this model are the following: 

0i,t+i^9{0i,t,ri,t) 

N 

r,,t+^ = l + J2Mrj,t-^) (37) 

where i e {1, • • • , N} denote the index of the i*'* rotator in the network. The function 
g{9, r) = 29 + b{r — 1) is the same as before and the real numbers Jij form a matrix J. The 
only requirement about this matrix is that its spectrum lies inside the unit circle such that 
it gives in fact a contraction of the radial variables onto = 1 for any i. On the other hand 
we can interpret J as the connectivity matrix of the network, meaning that the rotator j 
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influences rotator i iff 0. Moreover tliis interaction is signed, witli tlie interpretation 

tliat positive means activation and inliibition is described by negative J^. 

The model (37) is a straightforward generahsation of system (28), so that the explicit time 
evolution for (^j rj^t) can also be calculated. Furthermore, if J is diagonalisable, with 
P~"^ J P = dia^'je"''!"'"*'^!, • • • , e"'"^"'"*'^^} , then the stable complex susceptibilities (with 
respect to a perturbation of amplitude eX{6) in the Vj direction as considered in the single 
rotator model) can be expressed as follows: 



) = Y- 



1 1 



k 



So, it is readily seen that there are stable poles situated in the lower half-complex plane 
va. uj — ujk — ifJ^k and in uj = ujk — i(21og2 + /i^), each time e~^'''^'^^'^ is an eigenvalue of 
matrix J. Therefore, although in this model the unstable resonances are located only on the 
imaginary axis, depending on the connectivity matrix there can be an arbitrary number of 
stable resonances with non-zero frequencies uj^- The latter will not be visible in the power 
spectrum of the chaotic dynamics, but eqs. (38) shows that a periodic forcing of unit j 
at frequency U}. can produce a resonance in the linear response of unit i. As discussed in 
ref. [11], this property can be used to transmit a signal from j to i (possibly with amplitude 
modulation). Moreover we see that a necessary condition to make this transmision possible 
is that PikPkj 7^ 0. This condition is satisfied if the k*^ eigenmode of J has a nonzero 
component i in the canonical basis and reciprocaly if the j*'* vectors of that canonical basis 
has a nonzero component k when decomposed in the eigenvector basis. (Roughly speaking, 
an excitation of the node j excites the eigenmode /c, which excites the node i). This opens 
the possibility, for appropriately chosen connectivity matrices J, to get resonances which 
are specific to the pair ij. Note in particular that if the network as some modular structure 
such that the matrix P has a block diagonal structure then, obviously, there are resonances 
specific to each block. This structure is revealed by our excitation-response procedure even 
when one has no access to the explicit form of the dynamics. 



4 -3 Relaxation time longer than the mixing time 



To summarise, in the previous sections we have seen that the linear response between two 
variables (xj, Xj) of a uniformly hyperbolic system can be characterised by poles of its suscep- 
tibility Xij{ijj). A given pole uj — ujk — il-ik of this function is associated with two time scales, 
namely its resonance frequency cuk and its relaxation time /x^^. The latter gives an estimation 
of the decay time of oscillations Uk when the system is submitted to an impulse perturba- 
tion via its Xj variable. Moreover, in view of the decomposition Xiji'^) = Xij'*('-^) + Xifi'^) 
introduced in eq. (27), there can be two types of poles. The first class of poles are the ones 
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of x\j{u;). They coincide with the Ruelle-PoUicott resonances, which are also the poles of 

the power spectrum of a generic variable of the chaotic system. The knowledge of the closest 
pole to the real axis in this class provides the mixing time r^, defined as the inverse of 
its imaginary part. This is also the typical decay time of the correlation functions which 
can be interpreted as the "decorrelation" time of state variables of the chaotic attractor. 
In the basic examples treated above = l/log2. The second class of poles is formed by 
the poles of the stable susceptibility x\j\u). In particular their imaginary parts characterise 
the average decay time of an impulse perturbation along the stable direction, and thereby 
the average time to come back to "equilibrium" after this perturbation. For instance in the 
example of the previous Section, eq. (38) shows that the return time of < >t to equi- 
librium after an impulse perturbation in the stable direction rj is dominated by the largest 
Tij — maxfe(21og2 + {k chosen in {1,---,N such that Pi^P^^ 7^ 0}). Thus in this 

example r^ixing > Tij for all pairs ij and this corroborates the common belief that the return 
time to equilibrium of a physical system is given by its mixing time. 

This property is not true in general, however. The goal of the present Section is to support this 
statement by providing a generalisation of our simple model (28) where indeed the decay of 
response functions can be much longer than the mixing time. Moreover since this property 
was already observed in our previous study of neural network dynamics [11], the present 
example indicates also a possible mechanism to understand this non-classical behaviour. 

Let us consider the following system, which is an extension of system (28): 



^t+i = 2et + h{rt) {n - 1) + {h{rt) - 1) sin 9t 

ri+i = l + e-^(r-l) (39) 

with h{r) = e^^'^'^^^^^ and 77 > 0. Thus in the neighborhood of |r — 1| <^ i], the new term 
ip{rt) — 1) sin 9t is negligible and the system then behaves exactly as the basic model of 
Section 4.1. This is no longer the case when r — 1 increases, entaihng that h{r) decreases 
to 0. Figure 4 shows the 9 dynamics in the two limits r — 1 and r — 00. This model is 
hardly tractable analytically because it is hard to compute {9t,rt) as a function of (^qj^o)- 
But the dynamics can still be analysed by making use of the numerical method described in 
Section 3.2. Figure 5A shows a superposition of the susceptibilities |x6)r('^)| obtained from 
two numerical computations performed with two different values of rj. For 77 = 1 in eqs. (39), 
the numerical results can be well fitted with the analytical curve computed in the case 
r] — 00 (i.e. the Fourier transform of eq. (38) or, equivalently eq. (34) in the case of only 
one unit). Then the width of the resonance in u; = is dominated by log 2 corresponding to 
the inverse of the mixing time. On the other hand the narrower curve (with symbols "x" ) 
in fig. 5 shows an example of |x6ir('-o')| computed for small 77. In this case it is seen that the 
width of the resonance in a; = has substantially decreased. As a matter of fact it can be 
arbitrarily decreased by diminishing further rj. In the same time the height of the resonance 
peak increases. (This is not visible on the figure because their maxima have been scaled to 
the same height). Thus the susceptibility X6»r('^) diverges with 77 — > because one of its poles 
pinches the real axis in this limit. The reason for this behavior can be understood at least 
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Fig. 4. First return map representation of the 6 dynamics defined by eq. (39) in the two hmiting 

;( 




Fig. 5. (A) Modulae of the susceptibiUty Xori^) computed for the model (39) for different values 
of r]. The dashed line curve is the analytical prediction of |xer('^)| in the case rj = oo. The case 
r] = 0.05 (symbol 'o' ) is numerically computed with estimate (21) (e = 0.01, T = 2^^ and yet 
averaged over 200 samples). It shows a remarkable narrowing of the peak. (B) Impulse response 
functions X9r{t) computed from the inverse Fourier transform of X0ri^^)- In both cases (A) and (B) 
the height of curves r] = oo have been scale to the same amplitude than ij = 0.05. 

qualitatively as follows. As r — 1 increases from t], the system (39) looses indeed its uniform 
hyperbolicity. This is clear on fig. 4, where it is seen that it r — 1 ^ rj then the first return 
map for 6t acquires a slope 1 in abscissa ^ = 0. This means that a perturbation along the 
radial direction can induce a drastic slowing down of the "chaotic" rotation. The latter is 
resumed only when the radial perturbation has sufficiently relaxed. Thus in this case this is 
the decay time along the contracting direction which controls the return time to equilibrium. 
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Now one interesting physical consequence of this phenomenon appears by reconstructing the 
impulse response function Xer{t) from the inverse Fourier transform of Xer{^)- Figure 5B 
shows the comparison between the cases rj large (cf. eq. (34)) and rj small. In the first case 
the response of 9 decays as expected like 2~* which is the typical decay time of the correlation 
functions of this system. In the second case it is seen that the system responds with an impulse 
of large amplitude which moreover decays on a time of order 10 times longer than the mixing 
time. Therefore this example illustrates a novel feature of Ruelle's theory, compared with 
the common intuition of nonequilibrium physical (conservative) systems: when a dissipative 
system is pushed out of equilibrium, its return time to equilibrium can be substantially 
longer than the mixing time which could be predicted by looking at its correlation functions. 
This property merely reflects the possible violation of the fluctuation-response theorem in 
dissipative systems. 

The numerical simulations reported in this section concerned our single rotator which is a 
2 degrees of freedom system. Nevertheless we anticipate that the same type of phenomenon 
will show up in the network of rotators described by eqs. (37). In fact this behaviour of 
having a longer return time to equilibrium than the mixing time was also observed in the 
numerical simulations of our neural network dynamics [1 1] . 



4-4 Decomposing the linear response in general coordinates 

The model (28) is sufficient to show the existence of two types of responses in hyperbohc 
systems. However, the choice of variables (9, r) made to perturb this system and to observe 
its response is quite particular because its coincides with the local unstable and stable di- 
rections {eo,er}. So the decomposition of the response functions, or of the susceptibilities, 
into their stable and unstable parts is trivial for this choice. In this Section, we analyse the 
same model, but with observables which are the cartesian coordinates in the plane. The un- 
perturbed dynamics is thus the same but we change the way to perturb it and the variables 
to look at. The problem becomes slightly more complicated because each element of the the 
susceptibility matrix have both components, stable and unstable. The identification of these 
two contributions in a response function can scarcely be done explicitly in a general system 
because the projectors onto the local stable or unstable manifolds are not easy to compute 
for an arbitrary system. Nevertheless the present model is sufficiently simple to perform this 
decomposition analytically and to discuss some of its consequences on the resonances of the 
susceptibility. 

Let us consider the model (35) introduced before but now written in cartesian coordinates 
{x, y) of the plane: 

xt+i = F^{xt, yt, xt-i,yt-i) + ^^{t + 1, xt, yt) 

Vt+i^ Fy{xt,yuXt-i,yt-i) + eyCy{t + l,xt,yt) (40) 
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with the functions F = {F^, Fy) defined by: 



Fx{xt, yt, xt-i,yt-i) = R{\jx1 + y1, ^x1_^+yl_i) cos(^[arg(a;, y), ^ xl + yl]) 



Fy{xt, yt, xt-i, yt-i) = Ri^/xf+y^, ^Jxj_i+y^_i) sm{g[arg{x, y), V^fTyf]) 



(41) 



and the functions R and g have aheady been defined after eq. (35). The notation arg(a;,|/) 
means arctan(y/a;) for x > and arctan(y/x) + tt for x < 0. Thus in the absence of per- 
turbation (e^ = ey = 0), the standard change of variable x — rcos9, y — rsvaO gives the 
unperturbed model (35). 

Assuming without loss of generality that the mean values of the variables are zero in the 
unperturbed situation, < x >— =< y >, we are interested in the mean values < x >t and 
< y >t for small e and a prescribed form of the perturbation function ^(x,y). As already 
discussed (see footnote 5), in order to simplify the formalism we choose a perturbation the 
form ^ = XoF, with X = X^e^ or X = XySy. The goal is now to compute the susceptibility 
matrix Xjj(a;), or equivalently the corresponding response matrix Xijit) defined in eq. (25), 
as well as its decomposition in a sum of an unstable and of a stable part. To this purpose, 
by construction in the present model the perturbation XjBj {j — x or y) can be decomposed 
along its stable and its unstable projections, for example : 

XySy = Xy Slli 9 G,. + Xy COS 9 Sff (42) 

where is the radial unit vector (local stable direction) and is the tangential unit vector 
to the circle (local unstable direction). Thus from eqs. (25)-(27), and by using the equations 
VFj.e^ = drFi and VF^.ee — ^deF^ {i — x or y) , one can write the decomposition of Xxyi^) 
into the stable and unstable susceptibilities as follows: 



for t > (and Xxy = otherwise). The other matrix elements Xxxit) or Xyyi't) have an 
anologous expression. Let us recall that the invariant measure to be used in these integrals 
takes the following form (in polar coordinates) : ppidx) — rd9S{r — l)S{g — l)drdg with 5 
denoting here the Dirac distribution. This enables one to get the following results: 



Xxy{t) = xi'>{t)+X%Kt) 

/f COS 

PF{dx) sin 9 a,F^(x) A:^(x) + J p^idx) ——deF'^{x) Xy{x) 



(43) 



xrj(t) = < cos^t; d0[sm9oX,{9o)] > 
X^i{t)=drR' < cos 9t; cos^oX,(^^o) > 
X^xjit) = - < cos9t; de[cos9oX,{9o)] > 
xlc^{t) = -drg' < sm9u sm9oXy{9o) > 



(44) 
(45) 
(46) 
(47) 
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(We do not write the corresponding functions X^Ji^) X^Ji't) which are quite similar). We 
have only considered a ^^-dependence of the perturbation function Xi{6) (because a r depen- 
dence would be trivial due to the Dirac distributions). As it is expected the unstable response 
functions (44) and (46) are pure correlation functions, whereas the stable responses (45) and 
(47) are correlation functions multiplied by a factor related to the contracting dynamics. As 
already discussed this latter term is responsible for creating poles in the susceptibility which 
are not observable in the powerspectrum. 

Now, another point that we want to illustrate in this section is the following: as mentionned 
above, the decomposition xifit) + Xi]\'t) can scarcely be done in practise. In this case the 
resonances which can be observed arc induced from the maxima of the modulus |Xji(^)l- 
But the latter can be misleading to detect actual poles of the stable or of the unstable 
susceptibilty. To make this point clearer, we treat a concrete example. 

Let us consider again the perturbation function X{9) — 9{2t^ — 9) to be used in eqs. (44)- 
(47) and let us compute explicitly Xxx{t) = X^J{t) + Xxyi^)- The first term is given by the 
correlation function (44) which can be written as: 

7 . ff + |7r3 ift = 

X^:Jit) -jd9 cos(2*^)a,[e(27r - 9) sin 9] = l'^ ' ^ (48) 

27r ( (2t+ij2 + (2t_i)2 + (4t-i) j if ^ > 



and, as usual, xi"H^) = for t < 0. Next we compute the stable susceptibility Xxli^) 
evaluating the two factors of (45). The first factor S^-R* is worked out by using (36) and one 
gets 

sm ujq 



with the second factor being the correlation function: 



27r 



C{t) ^ Jd9 cos(2*^) cos(^) ^(1 -9)^ { 





. -27r ((2*^ + (2*31)2) ift>0 



(50) 



Finally the complex susceptibility is obtained by summing the stable and the unstable con- 
tributions. For example Xxx{t) — X^Ji^) + XxK'^) '^^^^ eqs. (48)- (50). 

In this example the poles of the complex susceptibility can thus be controlled by changing 
the parameters n and cuq which characterise the radial (thus here "transverse") dynamics. 
Figure 6 shows the modulus of the complex susceptibility |Xa;a;(^^)|) as well as the respective 
contributions of the stable and the unstable susceptibilities, for a given choice of parameters 
a = 0.8, a;o = 1. For this choice one notices that the resonances exhibited by these curves 
are all situated in a; = tt. So in this case the information given by the Fourier transform of 
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Fig. 6. Complex susceptibility x^xi^) computed from eqs. (48)-(50) with parameters 
a = 0.8, Wo = 1- (a) Modulus of the complex susceptibility, (b) Modulus of the stable (solid 
line) and of the unstable (dashed-line) susceptibility. 

a correlation function, would suffice to characterise, at least qualitatively, the susceptibility 
of the system. In particular, the maximum response of the system to periodic forcing occurs 
at the same frequency that the resonance of correlation functions. This situation changes by 
tuning parameter uq, however. Figure 7(a) shows a case where the resonance of the complex 
susceptibility is no longer situated in = vr, but near 7r/2. 

As already discussed, this "anomalous" behaviour cannot be understood in the framework 
of the fluctuation-response theorem, but well in the Ruelle's response theory of dissipative 
chaotic systems. On the other hand it can be remarked that the resonance of |x^a,(a;)| near 
7r/2 does not correspond to a stable nor to an unstable pole. In fact it results from the 
interference of two poles, stable and unstable, which appear as distinct resonances, which 
can be seen on the graph of 1x1^2(1^)1 and of Ixixli^)! represented on Fig. 7(b). 

So this simple example illustrates that the linear response of the system depends on the 
behaviour or the "stable" response function. In particular this can drastically differ from all 
those of correlation functions. Moreover, one sees that the position of the stable poles can 
be tuned by controlling appropriate parameters ruling the dynamical systems. 



Amplified echo of the impulse response 

We end this Section by reporting another property of the stable response which is an applica- 
tion of the analytical results (44)- (47). The following result is also linked to the interpretation 
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Fig. 7. Complex susceptibility x^xi^) computed from eqs. (48)-(50) with parameters 
a = 0.8, uo = 3. (a) Modulus of the complex susceptibility, (b) Modulus of the stable (solid 
line) and of the unstable (dashed-line) susceptibility. 



by mean function B(t) [eq. (33)] of the transient growth of response functions, as discussed 
in Section 4.1. 

The response function is a product of an exponential function and of a correlation 

function. If the latter is of amplitude of order (1) at time t, the result gives an amplification 
of 2*. To illustrate this property consider a perturbation of the form: 



X{9) = cos(2P + 1)^ 



(51) 



Then is easy to work out the correlation function < sin 6t sin ^o-^j/(^o)] > appearing in (47). 
It is found that if p > 1, the stable response function is 



Xxy (^) 



b{2P 



4(2 - e- 



H,p 



(52) 



Notice that in this case the susceptibihty x'^}{(^) oc e*^'^ has no poles. However, it is seen in 



this example that an effect of the stable direction is to produce an exponentially amplified 
"echo" of the perturbation. 
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5 Conclusion 



In the recent theory developed by Ruelle [4] the standard concepts of hnear response, sus- 
ceptibihty and resonance are revised and extended due to the dynamical contraction of the 
whole phase space onto attractors. More technically, the contraction rate of phase space 
leads one to replace the conservative Liouville measure (which is the standard measure in 
classical statistical physics) by the notion of SRB measure. Among others, this change de- 
mands to go beyond the standard frawework of the "fluctuation-response" theorem, which 
means that one cannot any longer deduce the statistical behaviour of the system just by 
looking at its correlation functions. For example, the powerspectra of classical (conservative) 
physical systems provide in general all the informations about the resonances of this system. 
For the dissipative systems considered by Ruelle, this information is incomplete, however, as 
they may be new types of resonances -called the stable resonances because they are related 
to stable manifolds- which cannot be detected in the spectra of correlation functions but 
in the complex susceptibility of this system. In previous papers we have exhibited a model 
of dynamical network with dissipative chaotic dynamics where we succeeded to numerically 
compute the susceptibility and show indeed the presence of new resonances in the suscep- 
tibility which are absent in the Fourier transform of the correlation functions. In order to 
gain more intuitions on what bring about these "stable" resonances and their consequences, 
we have considered in the present article simpler models which could be treated analytically. 
Moreover the analytical results could be used to test further our numerical method and to 
compare it to a previous method which existed in the litterature. 

We have first examined the 1-dimensional model 9t+i — 26t (mod 27r) which presumably 
is the simplest chaotic system with Lyapunov exponent log 2. Despite its simplicity this 
model is instructive to consider, as it is a non classical physical systems (having a purely 
expanding and non invertible dynamics). The linear response of this system can be worked 
out analytically and the result is a correlation function. In this case it is easy to show that 
this response decays exponentially with a characteristic time equal to the mixing time. Also 
the latter is clearly associated with a pole of the susceptibility and also identified with one 
Ruelle-Pohcott resonance. 

Next we have studied a class of 2-dimensional models describing simple rotator models. The 
angular dynamics is the same as previously but there is also a radial dynamics which contracts 
the phase space on the unit circle (fig. 3). The resulting system has thus a non- fractal chaotic 
attractor, which moreover is uniformly hyperbolic. This elementary model could be used to 
illustrate several new properties associated with the dissipative dynamics. First it is easy to 
show in it the presence of new poles whose real parts correspond again to zero-frequency 
resonances but imaginary parts are associated with the relaxation time of the contracting 
dynamics. In fact it is seen that the cross-response X6>r(<^) has a pole which combines the two 
characteristic times (mixing and contracting) in a linear combination of both. This property is 
demonstrated to exist in a class of more general and higher dimensional systems, as presented 
in the Appendix. In the rotator model one could also describe the stable resonance of Xeri'^) 
by means of the error function B{t), eq. (33), showing how the analytic properties of the 
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perturbation function X{6) regulates the amplitude growth and decay of the response. Next 
we extended the calculations obtained for a single rotator to resuh.s for a network of coupled 
rotators. This generalisation pointed out to the possibility to get stable resonances wich are 
specific to the pair ij in the network and having arbitrary non zero frequencies. Both of these 
features can be shown to be intimately linked to the spectral properties of the connectivity 
matrix of the network. This properties would be worth to analyse further as they are also 
related to the ability of the system to propagate signals between units of the network [11]. 

We have also considered an extension of the rotator model in another direction which allows 
one to discuss the relation between the return time to equilibrium of the system after an 
impulse perturbation and the mixing time. These characteristic times coincide when there 
is no attractor in the phase space but here we proposed a mechanism to get a relaxation 
time arbitrary long in case of transverse perturbation to the attractor. Roughly speaking 
it is enough that the transverse dynamics becomes less contracting with the distance to 
the attractor leading possibly to the loss of hyperbohcity for critical perturbations. More 
precisely we have implemented such mechanism in model (39) which presumably is beyond 
the scope of analytical treatment, but could be dealt with numerical simulations. These 
numerical results indeed showed the narrowing of stable resonances following our qualitative 
predicitions. 

Finally we have analysed the simple rotator model in cartesian coordinates, illustrating a rare 
case where the decomposition of the susceptibility Xiji^) = Xifi^) ~^Xij\^) respectively 
stable and unstable components can be performed analytically. This example showed in 
particular that in some cases a maximum in the amplitude |Xij(i-i-')| may not be interpretable 
simply as a stable or an unstable pole, but rather as an interference effect of both. 

The results obtained in this paper gives more insights in our previous work concerning the 
linear response of a chaotic network model [8]. This preceding study had allowed us to nu- 
merically demonstrate the existence of the new resonances predicted by Ruelle in a model of 
neural network. Thanks to a numerical method devised to compute the susceptibility and to 
extract also some poles of it, these new resonances were observed with the unexpected fea- 
tures to be often narrow and pair ij-dependent. This last property enabled us also to propose 
a method of signal transmission bewteen the nework units by means of amplitude modulation 
of harmonic signals locked on the resonance frequencies [11]. Now the present work permits 
us to understand some mechanisms able to create such resonances with arbitrary frequencies 
and arbitrary widths, and the possibility to get pair-selective resonances. In future works we 
will examine further how the control this pair selectivity and how pole positionning can be 
achieved so that we can propose more advanced schemes for communication between units 
within a chaotic network. We anticipate that such schemes would be quite useful not only in 
the field of neural networks but also for applications in other biological networks (e.g. gene 
regulatory networks) or in communication networks. In particular we will examine also ways 
of transmitting trains of periodic pulses, which are prevalent in the considered applications. 
On the other hand, to demonstrate the practical relevance of these new ideas we are working 
on experimental devices with dissipative chaotic dynamics in order to produce experimental 
evidences of the present linear response theory [17]. 
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From the theoretical side it would be worth to deepen the consequences of the fractal (cantor- 
like) structure of generic strange attractors, a feature which was not considered in this paper. 
A generalisation of our rotator model would be to replace the contracting dynamics along the 
radial direction by a baker map. To this purpose a good model candidate to be investigated 
in the future could be the nonlinear multibaker map analysed in [18]. This reference considers 
a class of perturbed multibaker map mimicking a field driven, thermostated random walk 
for which a non trivial SRB measure can be analytically computed to first order in the field 
parameter. This system exhibits indeed a fractal structure in the stable direction. Moreover in 
some limit, this system becomes nonhyperbolic, in a similar way to the example numerically 
treated in Section 4.3. 

Finally, many questions addressed above would be worth to be studied beyond the hnear 
response regime, which is in the scope of our present numerical techniques. 



A Some generalisation in n-dimensional space 

In this appendix we analyse further the way the impulse response function Xiji't), defined 
by eq. (25), can be decomposed as the sum Xijit) = Xifi^) + Xi]\'t) [sq- (27)] in a higher 
dimensional space. The unstable contribution xtj\t) this sum has been dealt with in 
details by Ruelle who introduces the concept of unstable divergence. In what follows we 
recall this concept and we present an attempt of dealing with the stable part Xi]\t)- To 
this purpose we introduce hypothesis on the underlying dynamical system so that we can 
generalise some conclusions obtained in this paper with the simple rotator model. 

The unstable part Xij{t) has been studied by Ruelle in [4] in the case of a general observable 
A. When A = Xi the result of this analysis can be expressed by a correlation function which 
can be written as: 



In this equation the differential operator div" is called the unstable divergence, i.e. the usual 
divergencJ^ operator restricted to differentiating the field X" only along the unstable di- 
rections and computed with respect to the invariant measure pp{d'x) (which is absolutely 
continous in these directions). We propose to retrieve the result (A.l) in the particular case 
where we can assume a change of variables 



^ The divergence div^^X = ^div(CX) of vector field X computed with respect to the invariant 
measure p{d-x) = C(x)(ix can be defined as the Lie derivative Lx{p) = (div^X) p. 



X^\t) = -j p,{d^) i^*(x) div"X«(x) 
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such that u and s are coordinates representing respectively the unstable and the stable 
manifolds of system Xt_|_i = F(xt). Moreover let us suppose that coordinates (u, s) can be 
characterised by a diagonal metrics: 

m n 

||rfx|p = ^5(fc(u, s)(iMfc + ^/ife(u,s)(i4 {m + n^N) (A.3) 

fe=i fe=i 



with = Wq^W^ and = II^^IP- These hypothesis are a direct generalisation of the 
rotator model studied in Section 4. In that case the coordinates (u, s) are simply the polar 
coordinates {0,r). 

Using these coordinates, we assume also that the invariant measure pF{d:x.) can be factorised 

as 

PF{dx) = p^''\du) p^'\ds) (A.4) 



where the measure in the stable directions, p^^\ds), is typically singular, and p*^") is absolutely 
continuous along the unstable directions with a density C(u). So the invariant measure pF{dx) 
will be also written as: 

PF(dx) = C(u) \fgin^dui ■ ■ ■ dumP^^\ds) (A.5) 



where ^ is the determinant of the metric tensor (restricted to the unstable coordinates). 
Now we compute the unstable contribution of the susceptibility : 



it) = / pF{d^)VF^{^) . X(.")(x) (A.6) 



where X^""* is the projection on the unstable space of vector (the perturbation along the 
coordinate Xj). We can assume that X^"^ decomposes as: 

Xf (u,s) = ;^X,,(u,s)-^ (A.7) 
fe=i ^ 



Then by noticing that VF/(x) • X^"^ = SfcLi ^-^jk the unstable susceptibility can be 
written as: 

(t) = J P^'^ ids) CVgdUl--- dUrr. 9^^.-fcJ 



Integration by part on variable Uk simplifies the equation because the integrated term can 
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be assumed to vanish. Therefore one obtains: 

X^\t)^-Jp^'\ds)dm---dumF!f:^[C,/9X,,] (A.8) 

k=l ^ 

Finally the equation (A.l) is retrieved by defining 

When C = 1 the righthand side of this equation is the standard expression of the divergence 
expressed in curvilinear coordinates. The upperscript u indicates however that the differen- 
tiation acts only with respect to the coordinates u of the unstable manifold (see footnote 
7). 



The stable response function 

We focus now on the stable contribution of the decomposition (42) of the susceptibihty. This 
term is thus given by: 

xff W = / PF{d^)VFl{^) . X«(x) (A.IO) 



where X^-*^ is the projection on the stable space of the perturbation X^ . Let us remark that 
X^"*^ can be computed in terms of the function ^ as: 

X?^ = ^.E7^I^# (A.11) 



So the integral for xij becomes sightly more exphcit: 

.g'w^/p<"'(..)/"v.)|:(i£A-g 



The reason why we cannot transform this expression into a correlation function, as was done 
for the unstable part, is that here we are not allowed to integrate by part on variables s^. 
Indeed the invariant measure is typically singular in s (there is no density function in p^^^). 
At this stage one cannot proceed further in the development of xfjit) unless additional 
hypothesis are made on the dynamics. 
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In what follows, simplifying hypothesis will be added on the dynamics in order to factorise 
the integral (A. 12) into separate contributions of respectively the stable and the unstable 
coordinates. Although these assumptions are compelling, they are satisfied in the case of our 
simple rotator model, thus generalising the properties obtained in this framework. 

So let us consider four factorisation hypothesis as follows. The first one is that the dynamics 
Xt+i = F(xt) decouples in terms of coordinates (u, s) and becomes 



ut+i = U(ut) 

St+i = S(sO (A.13) 

with some functions U and S. The second assumption is that the change of variables from 
(u, s) to X can be written as the product 

Xj = ^j{u, s) = cpj{u) 7jj{s). (A. 14) 



By combining these two assumptions the t-th iterates of dynamics can factorise as i^/(x) = 
93j(U*(u)) ?7i(S*(s)). The third hypothesis is that the metrics functions hk defined in eq. (A. 3) 
can be also factorised: 

hk{n,s) = ak{s)(3k{n) 



Then, assuming finally the following choice for the perturbing function Xj(x) — r(s)'0(u), 
the stable susceptibility can be factorised as: 



xifit) = t [j '^^^^(^^)^|r(S*(s))g^(s)r(s)) (/ pW(d«)^^,(U*(u))vP,(u)^(u)' 



akdsk dsk J \J pk 

Y.Z,^k{t)C,,k{t) (A.15) 



k=l 

n 



k=l 

namely 



In this last expression of xfj (t), the Cijk{t) are correlation functions of the unstable dynamics. 



Ci^k{t) =< ^i{nt)-, <^,-(uo)|^ > (A.16) 



whereas the Zijk{t) are not correlation functions because it cannot be expressed as (14). 

In order to simplify the notation, now we consider the case where there is only one stable 
direction {n — 1). Then eq. (A.15) reduces to: 

X^{t) = Z,,{t)Q,{t) (A.17) 
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For hyperbolic mixing systems the correlation function Cij{t) can be decomposed as a sum: 

oo 

k=i 

whose Fourier transform is readily deduced as: 
CiAu) = .^'-^^ , h c.c. 

k=l 

This function possesses complex poles for u; — ±0;^ — i/ik- These are the Ruelle-Pohcott 
resonances already mentioned and they do not depend on the choice of the pair ij. For 
hyperbolic mixing systems one gets < fii < ^ ■ ■ ■ and fii^ is the mixing time (cf. 
Section 4.3 ). Thus the Fourier transform of (A. 17) becomes the convolution product: 

xffH = ^ / ^^Jf^^i^') (A.18) 

where Zij{u) is the Fourier transform of Zij{t). The latter is not a correlation function 
but it can be assumed to exponentially converge towards since it is related to the stable 
dynamics. So, due to the singularities of Zij{uj) the complex susceptibility x,[^^(cj) can acquire 
new poles distinct from the Ruelle-Policott resonances, and the latter can depend on the pair 
ij as mentionned in Section 4.2. For example let us consider the simplest case for Zij{t) where 
it is a function of the form: 

Zij{t) = aij e-(*^+^P*^)* + c.c. > 0). (A. 19) 

Then the complex susceptibility takes the expression: 

<>{')(uj) = 'Y ( ^'^^ h ^^^^ \+cc fA20) 

Therefore, in this simple case where Zij{Lo) can be characterised by only two poles q-ij ±ipij, 
it is seen that the stable resonances arc obtained by splitting the Ruelle-Policott resonances 
in any combinations of the form icufcipj^ — i(/i,fc + gjj). Notice that the real parts of these new 
poles are bounded from below by ^i. So in this framework the time of return to equilibrium 
is the "mixing" time. 

Extending the example of eq. (A. 19), new resonances can be created for each poles of Zij and 
for each terms (i.e. each stable directions) in the more general expression (A. 15) of Xiji'^)- 
Therefore we conclude that the existence of these poles enables the susceptibility to possess 
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resonances which possibly are specific to the pair ij. This means that a periodic forcing of 
one degree of freedom j of the system can resuh in principle in a resonance in another degree 
of freedom i whose frequency is specific to it. This is precisely this possibility of ij-dependent 
resonances which is exploited in [11] in order to perform selective transmission of signals in 
a chaotic networks of interconnected units. 
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